Diffusion Monte Carlo: exponential scaling of computational cost for large systems 
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The computational cost of a Monte Carlo algorithm can only be meaningfully discussed when 
taking into account the magnitude of the resulting statistical error. Aiming for a fixed error per 
particle, we study the scaling behavior of the diffusion Monte Carlo method for large quantum 
systems. We identify the correlation within the population of walkers as the dominant scaling 
factor for large systems. While this factor is negligible for small and medium sized systems that 
are typically studied, it ultimately shows exponential scaling. The scaling factor can be estimated 
straightforwardly for each specific system and we find that is typically only becomes relevant for 
systems containing more than several hundred atoms. 



PACS numbers: 02.70.Ss, 71.15.Nc, 31.15.-p, 



I. INTRODUCTION 

Today's scientists can choose from a wide range of 
computational methods for the simulation of quantum 
mechanical systems. These range from highly efficient 
semi-empirial methods to density functional methods - 
offering a practical compromise of efficiency and accu- 
racy - all the way to very accurate quantum chemical 
methods. Besides these deterministic methods, various 
stochastic quantum Monte Carlo (QMC) methods are 
gaining ground, offering exact handling of many strongly 
correlated systems and scaling up to system sizes that 
are out of reach for the deterministic competitors. 

The two major arguments that are typically brought 
up in the advocacy of QMC are the excellent paralleliz- 
ability and the good scaling behavior. Depending on the 
QMC variant that is chosen, the collection of statistical 
data points can be performed in parallel with little to no 
communication, making the method well suited for high 
performance computers of all architectures. The scaling 
behavior depends greatly on the details of the system 
and the method, but it is generally found to be signif- 
icantly better than that of quantum chemical methods, 
and linear scaling algorithms have been reported^— 

A commonly used method for the ab initio sim- 
ulation of electronic structure is diffusion Monte 
Carlo (DMC)j£~— typically using the fixed node 
approximation^ For this method, the bulk of the com- 
putational effort is spent on the repeated evaluation of a 
trial wave function for electron positions that change step 
by step, one electron at a time. The trial wave function 
is usually expressed as a Slater determinant^ of single 
electron orbitals, multiplied by a Jastrow factor— to ex- 
press electron correlations. For single electron orbitals 
expressed as maximally localized Wannier functions^ 
the local energy can be reevaluated in constant time after 
a single electron move, leading to an O(N) algorithm for 
a complete time step of one configuration^— As a further 
refinement to this, trial wave functions for DMC calcu- 
lations are today commonly expressed in a blip basis^ 
which can be evaluated very efficiently. 

In contrast to deterministic methods, however, the 



computational cost of a Monte Carlo (MC) simulation is 
meaningless without specifying the statistical error that 
is achieved. Deterministic methods typically have sys- 
tematic errors that are either intrinsic or depend on pa- 
rameters that do not scale with the system size. The sta- 
tistical error of MC simulations on the other hand scales 
very simply with the inverse square root of the run time 
while the scaling with the system size is a nontrivial issue 
that depends on details of the method and the system of 
study. Though the unfavorable scaling of the statistical 
efficiency of DMC has been demonstrated before,— it has 
- to our knowledge - never been studied systematically. 

In this paper, we will present a systematic study of 
the scaling behavior of QMC calculations aiming for a 
fixed statistical error bar. The main focus will be on 
the DMC algorithm including branching and population 
control as described by Umrigar et al. jiS other variants 
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FIG. 1: (color online) Scaling of the effective population size 



= 1.5, see Eq. JT3J] 
The target popu- 



Np" P [ see Eq. ([8J] in a sample system [a 
with increasing number of atoms A^tom- 
lation size is fixed, the true population size N pop fluctuates 
around a slightly lower average (see text). The "error bars" of 
iVp p visualize the increasing population fluctuations crjv pop . 
The effective population drops exponentially, due to increas- 
ing correlations within the population. 
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will be briefly discussed as well. The statements that we 
will derive are expected to hold for DMC calculations in 
general, but to simplify understanding, we will consider 
a "typical" system made up of N similar constituents 
which we simply call "atoms" . This could be, for exam- 
ple, a crystal in a simulation cell made up of N primitive 
cells, a cluster of N atoms or a large organic molecule of 
N comparable groups. 

We will begin by deriving several general quantities 
and continue by demonstrating these in the case of a 
simple model of N independent hydrogen atoms. From 
this model we can numerically extract the missing pieces 
of the scaling behavior, allowing a quantitative estimate 
of the scaling limit for arbitrary systems. This limit will 
then be discussed for a number of different sample sys- 
tems. 



II. SCALING OF COMPUTATIONAL COST 

The total computational cost of a DMC calculation 
(optionally split over a number of parallel CPUs) is 
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where N stcp is the number of steps in imaginary time, 
7Vp p is the average population size and t step is the CPU 
time needed for one single all-electron step per configura- 
tion. Using a so-called "linear scaling" QMC algorithm^ 
each all-electron move scales as 

t stc p oc N, 

assuming that the evaluation of the Slater determinant 
dominates the computational cost. For the moment, we 
assume that population fluctuations are negligible and 
-/Vp p can be treated as an external parameter. The in- 
fluence of the population control will be discussed later 
on. 

The standard error of the total energy can be expressed 

as 
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with the constituents explained in the following: 

The raw DMC variance (J^mc 1S the variance of the 
local energy of individual configurations over the whole 
simulation. Being based on the mixed estimator, cr| mc 
may deviate from the variance a vmc of the trial wave 
function obtained in a VMC run. For typical systems, 
however, we find that <r| mc ps er^ mc . Assuming that the 
trial wave function of the whole system can be optimized 
to about the same quality as that of a single constituent 
N, this variance scales as 



°dmc K N , 



(3) 



since the local energy is dominated by the sum of N 
independent atomic local energies. 



The (integrated) correlation time of a scries of data 
points Xi is given byi^ 
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where (•), denotes the arithmetic mean over the index i. 
T corr in units of the time step r s t p describes the factor 
by which the number of steps N ste p S has to be divided 
to correct the error bar of the result for serial corre- 
lation. To obtain T corr /r s t e p, as an alternative to com- 
puting Eq. directly, one can also use the rcblocking 
method.— Though the system may have various corre- 
lation time scales, some of which depend on the system 
size, we find that the integrated correlation time that is 
responsible for the reduction in the resulting accuracy 
is dominated by the shortest correlation times which de- 
pend on local properties, such as the kind of nuclei in the 
system, but not on the size N. 

The population correlation factor Xpop ^ 1 captures 
the inefficiency of the process due to population corre- 
lation and fluctuation. We will treat this factor as an 
unknown quantity for the moment and discuss it in de- 
tail afterwards. 

Using Eqs. ([T]) and along with the discussed scaling 
laws, we can express the scaling of the total computa- 
tional cost as 
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so we see that - apart from the factor Xpop ~ DMC is in 
fact a constant scaling method if a fixed standard error 
per atom 5E at om = 5Etotai/N is required, as it is the 
case for example in the study of long-ranged correlations 
in periodic systems. Of course, memory limitations or 
implementation issues will limit the size of computable 
systems. Within these limitations, however, the constant 
scaling behavior is not so surprising, considering that for 
collecting statistical data, it docs not make any difference 
whether you simulate N weakly interacting systems in 
parallel or a single system N times as long. Both result 
in the same statistical error for the same computational 
cost. 



III. POPULATION CORRELATION 

For small enough systems, the factor x pop is close to 
one, which may be the reason why, to our knowledge, 
a systematical study has never been attempted before. 
When scaling up the system size, however, population 
correlation becomes important and we need a better un- 
derstanding of its scaling. 

The DMC algorithm is based on a drift-diffusion pro- 
cess with branching and killing of configurations due to 
fluctuations in the local energy. A freshly branched pair 
of configurations is identical and thereby fully correlated. 
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In the following drift-diffusion process, it takes some time 
to decorrelate, leading to a fluctuating amount of corre- 
lation within the population at any time. 

We consider a DMC run over iV s te P time steps i. We 
will first consider a simplified model with constant popu- 
lation of iVp p configurations p, each having a local energy 
E p . A generalization including population fluctuations 
will follow in the section below. 

An effective population size N p ^ p can be defined as the 
number of configurations that would result in the same 
variance of the average as the correlated population. For 
a long DMC run, the raw DMC variance can be estimated 
from the averages over all configurations at all time steps 
as 
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size Wi = N pop for each time step i. The resulting total 
energy average of a DMC run is then 



Etat = 



iV B tep 

£ 
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In more sophisticated schemes, Wi and N pop may be de- 
coupled. To estimate the variance of this weighted aver- 
age, we can split off the correlation time into a factor and 
use the estimator of the variance of a weighted average. 
Viewing the local energies E p as random variables and 
the weights Wi as constants given by a long DMC run, 
we can express this as 



var[£'tot] = 



wf var 



Tstep (Ei^) 2 i= l 



while the variance of the population average is defined as 

2 



a 



pop 



In the same interpretation of E p as random variables, we 
can replace the averages over time steps in Eqs. ([5]) and 
([7]) by statistical expectation values (•) and write for each 
single time step 



with the averages abbreviated as (•) = J2p=T '/^pop 

and (■). = EtT P V^tep. 

In the case of an uncorrelatcd population, we would 
find CTp 0p = o"dmc/-^pop' so we can measure the amount 
of correlation by defining an effective population size as 
the ratio 
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Substituting this into Eq. ([9]) results in a sum over ex- 
pectation values, so the var [E tot ] itself can be written as 
the expectation value of a single expression which can be 
expressed as a product 
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where N p ^ p ^ N pop with equality only in the case of a 
completely uncorrelated population. By collecting the 
necessary data during a DMC run, iVpo p can be com- 
puted at negligible cost. Fig. Q] displays the scaling of 
the effective population size with increasing system size 
for a sample system. 

Note that the average population N pop is typically 
slightly lower than the target population iVp^p 60 *, because 
the population control implemented in CASINO uses the 
linear average ^££) of the energy instead of the exponen- 
tial average In (exp {E 1 ^)) p which determines the actual 
growth of the population. Apart from reducing the pop- 
ulation size, this has no effect on the statistics or the 
result. 



IV. POPULATION FLUCTUATIONS 

Keeping the population size completely fixed as we had 
assumed in the previous model gives rise to a population 
control bias. To reduce this bias, the population size 
N pop is allowed to fluctuate and only weakly controlled. 
The average over all time steps then needs to be weighted. 
In the simple population control mechanism considered 
here^ the weights are simply defined by the population 
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with generalized expressions for the effective step number 
and population size 
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where the quadratic appearance of the weights makes the 
effective population size sensitive to population fluctua- 
tions as well. 

To estimate the statistical efficiency of the DMC algo- 
rithm, it is most useful to combine both quantities into 
the definition 



Xpop 



Nstcp x N pop 
Nf ep xN^ p 



(10) 



For a DMC run that is sufficiently long that the set of 
weights Wi is a good representation of the statistical dis- 
tribution, the quantity Xpop is an unbiased estimator. 
It corresponds exactly to the quantity in Eq. (j^Jl and 
remains directly proportional to the total CPU cost ac- 
cording to Eq. ([5])- 
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V. ASYMPTOTICS OF x P o 



VI. MODEL DMC PROCESS 



The asymptotic behavior of Xpop for weak population 
correlation can be derived by a few simple arguments. 
Assume, for a moment, the local energy E\ oc of configu- 
rations to be normally distributed as 
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A single configuration with E\ oc < will branch at a 
rate of —E\ oc . Integrated over the distribution of E'loc, 
this leads to a branching rate per configuration of 
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It is, of course, known that the true distribution of the 
local energy is far from normali^I However, we can signifi- 
cantly relax the previous assumption, considering that we 
essentially obtained the ratio between standard deviation 
and mean absolute deviation which holds approximately 
for a wide range of distributions.— 

Starting from an initially uncorrclated population of 
size iVpop, the population after branching is A pop + 1 
with two identical configurations. The population mean 
is equivalent to that over A pop — 1 correlations of single 
weight and one of double weight which has the variance 
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For Ap 0p 3> 1, the effective population size after the 
branching is therefore (A 7 ^ )' = A pop — 1. 

To keep the population stable, the branching and 
killing rates have to be equal. For a weakly correlated 
population, it is trivial to see that a killing event reduces 
the (effective) population by 1 as well. 

After branching, the two copies evolving independently 
take an effective time of t co „/2 to decorrelatc. (Since 
T corr measures the amount of statistical data that is lost 
due to serial correlation, which is exactly the quantity 
that we want to measure for population correlation as 
well.) 

For Tbranch Tcorr j the effective population is reduced 
by 1 at a rate of 2 A, 
T corr /2. On average this gives 



p°p r branch an d restored to Ap p within 
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or, since population fluctuations can be neglected, 
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To study the dependence of Xpop on the system pa- 
rameters beyond the perturbative regime, we have im- 
plemented the full DMC algorithm on top of a minimal 
model of a correlated diffusion process. Each configura- 
tion is reduced to a single random variable with a simple 
exponential autocorrelation so that ffdmc and r corr are 
free parameters. The value of the variable is used di- 
rectly as the local energy for the branching process. 

The population control mechanism described by Umri- 
gar et al. introduces an additional parameter r CO rer with 
the dimension of time (implemented in CASINO as pa- 
rameter cerefdmc, used in updating the reference en- 
ergy E rc i; r cere f/T s to P corresponds to g in the original 
publication^). To avoid frequent population instabili- 
ties for extreme parameter settings, we have restricted 
the population size to a window around the target pop- 
ulation and decoupled the total weight from the total 
population size outside of this window, allowing recovery 
from explosions or starvation without introducing any 
additional bias. All the results presented below are in 
the regime where this mechanism has no significant im- 
pact on the efficiency. 

With each run, one obtains the factor Xpop as a 
function of the parameters <7dmc, Tcorr, T CG ref and r s tep- 
The result must be dimensionless, reducing the num- 
ber of relevant parameters by one. Furthermore, 
one is interested in the limit r s t e p — > 0. We find 
that for T stC p < 0.1 x min(r corr ,r C eref, 1/cdmc), the inef- 
ficiency factor Xpop becomes practically independent of 
T s tep in all cases. We can combine the remaining free 
parameters into erdmc Tcorr and T coro f /r corr , leading to the 
results displayed in Fig. [2] 

Most significantly, we find that Eq. (fTTj) is not only 
confirmed in the perturbative limit, but its exponential 
continuation gives a strict lower limit for the inefficiency 
factor 



Xpop ^ exp |o"dmc T corr / ^/2,7T 



(12) 



where the deviation from this exponential depends on the 
ratio T cere f /r corr . 



VII. HYDROGEN SAMPLE SYSTEM 

To demonstrate our result in a real calculation, we 
have performed various DMC runs using the CASINO 
program.— We chose a system of A hydrogen atoms 
placed several thousand atomic units apart to make them 
effectively independent. As a trial wave function, we used 
the exact ground state with a detuning parameter a and 
an additional term to satisfy the Kato cusp condition^ 
centered on each hydrogen atom 
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FIG. 2: (color online) Inefficiency factor Xp°p computed for 
a model DMC process (see text) in dependence of the two 
relevant parameters crdmci~corr and r C eref /Tcorr ■ The exponen- 
tial law extrapolated from the perturbative limit is found to 
give a strict lower limit for Xpop- Solid circles refer to stan- 
dard DMC with branching. Hollow circles refer to minimal 
stochastic reconfiguration MC (SRMC).— ^ 



We performed a large variety of runs on this model sys- 
tem with system sizes N <G {1,2, 4, 8,..., 64}, and detun- 
ing parameters a e {1.1, . . . , 3.0} and target population 
iV pop = 200. The DMC time step was set to r stcp = 0.02 
in all cases. 



The variance cr^ mc 



was found to be equal to a^ mc 
within the statistical error in all cases. In each case, the 
population correlation factor Xpop was determined from 
the variances using Eq. (fTUj) . 

Obtaining a precise value for the correlation time r corr 
takes an extremely large amount of data in either of the 
two methods described above. For a reasonably precise 
value, we performed a very long DMC run (-/V s to P > 10 7 ) 
on a single atom for each value of a. Several tests on 
larger systems confirmed that the same value holds for 
larger numbers of atoms N. Since r corr is independent 
of T s top only if r corr ^> T step , we performed these runs for 
the same step size as the main calculations. 

As in Fig. [2] we plot the population correlation factor 
Xpop from all our calculations as a function of the prod- 
uct OdmcT'corr and again find the exponential lower bound 
described by Eq. (|T2l , as displayed in Fig. [3] 
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FIG. 3: (color online) Data from many different calculations 
on a model system of TV independent hydrogen atoms with a 



detunable trial wave function [Eq. (|13|l ]. <7dmc an d Xpop were 
directly obtained from each run. r CO rr was determined from a 
single, very long run for each type of atomic wave function. 
Population control is kept at the default f CC rcf = 1. 
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or worse. So, even if population correlation may not be 
an issue yet for most applications, it will eventually lead 
to an exponential scaling of the cost. The factor X can 
be reduced by optimizing the wave function, but the gain 
that is possible with reasonable effort is very limited. 

Table fl] lists a selection of sample systems showing the 
size at which the exponential scaling becomes observable. 
The integrated correlation time r corr [via Eq. (j4])] and the 
raw variance cr\ mc were computed for very small systems 
and Eq. (jT3J) was then used to estimate the size at which 
a comparable system would show significant population 
correlation. All values should be understood as rough es- 
timates. The trial wave functions were either taken from 
a library of examples or optimized with moderate effort. 
Further optimizations could certainly reduce <J^ mc and 
thereby shift the onset of significant population corre- 
lation. Typically, however, significant effort is necessary 
even for minor improvements using optimizations beyond 
the standard Jastrow terms. 



IX. ALTERNATIVE VARIANTS OF QMC 



VIII. ANALYSIS OF VARIOUS SAMPLE 
SYSTEMS 

The exponential law in Eq. (|12[) has severe implications 
for the scaling of the DMC method. Following Eqs. §3§ 



To this point the discussion was centered on the con- 
ventional DMC algorithm including drift and branching. 
In the following, we will briefly discuss a number of al- 
ternative QMC algorithms in view of the population cor- 
relation scaling. 

First, it is clear that population correlation can only be 
caused by some form or branching. The variational MC 
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atoms (ae) 




Tcorr 


2 / 

o- dmc /atm 


Xpop — 2 


He 




0.5 


0.0044 


2700 atoms 


c 




0.4 


0.16 


140 atoms 


Ar 




0.04 


8.0 


250 atoms 


molecules (ae) 




°dmc/ mlc 




H 2 




0.1 


0.58 


550 molec. 


CH4 




0.3 


0.24 


120 molec. 


C2H4 




0.4 


0.51 


38 molec. 


SO2 




0.06 


7.5 


105 molec. 


crystals 






Odmc/ atom 




diamond (pp) 




0.15 


0.23 


630 atoms 


diamond (ae) 




0.1 


2.3 


133 atoms 


graphite (pp) 




0.3 


0.20 


135 atoms 


silicon (pp) 




0.4 


0.052 


328 atoms 


electron gas 






2 / i 

o-dmc/ elec. 




3d crystal (r„ 


= 1) 


0.2 


0.26 


193 elec. 


3d fluid (r s = 


5) 


5 


4.2 x 10" 4 


330 elec. 


3d fluid (r s = 


10) 


16 


5.1 x 10~ 5 


242 elec. 


2d crystal (r a 


= 1) 


0.4 


0.038 


570 elec. 


2d fluid (r„ = 


1) 


0.3 


0.033 


1154 elec. 



TABLE I: Estimated values for various sample systems. The 
last column gives the system size based on Eq. ()12|) at which 
the population correlation becomes significant with Xpop — 
2. Beyond this size, the DMC method must be expected to 
become exponentially inefficient. The first three categories 
are based on either all-electron (ae) or pseudopotential (pp) 
wave functions with optimized Jastrow factors. All numbers 
should be understood as rough estimates based on moderately 
optimized trial wave functions. Reducing Od mc by further 
optimization will shift the onset of the inefficiency by the 
same factor. 



(VMC) algorithm, which samples an explicitly known 
wave function, clearly does not have this feature. A 
variant of DMC with branching switched off (sometimes 
referred to as "pure" DMC) also features a completely 
uncorrelated population. If all configurations are fixed 
to the same statistical weight, this process produces the 
same distribution of configurations and thereby the same 
total energy as VMC and can therefore be seen as a vari- 
ant of the former. 

If, on the other hand, each configuration in a pure 
DMC run carries a statistical weight evolving with the 
fluctuations in the local energy, the effective population 
size is reduced in the same way as it would be 

when branching were allowed, with the only difference 
that decorrelation does not happen and the method be- 
comes exponentially unstable with simulation time.— 

A number of variants of the DMC algorithms keep the 
population size fixed and include branching in form of 
stochastic reconfiguration, duplicating some configura- 
tions and deleting othersi 23 ' 24 i 27 While the population 
is fixed, the total weight is allowed to fluctuate indepen- 
dently and the population control is replaced by a weight 
control mechanism. 

Our definition of the population correlation factor x pop 
in Eq. (jTT)]) is already kept general enough to capture the 
effects of weight fluctuations within the population along 
the correlations within the population and to capture the 



fluctuations of the total weight along with the population 
fluctuations. Tests on several variants of the branching 
strategy confirmed that these have no influence on the 
exponential lower bound of Xpop but only affect how far 
the actually measured Xpop exceeds the predicted expo- 
nential scaling (sec Fig. [5]). 

One remaining option to limit Xpop is the use of strong 
population control on a small population, accepting a sig- 
nificant population bias. The extreme case of this strat- 
egy would be a single walker with the weight renormal- 
ized after every step, leading exactly to the VMC distri- 
bution. We can, therefore, tune between exponentially 
scaling statistical inefficiency and a population bias that 
ultimately leads to recovering the VMC algorithm, which 
- as we know - does not suffer from exponential scaling. 

Variants of QMC such as path integral MC (PIMC) 2 ^ 
or reptation Monte Carlo (RMC)^ are somewhat re- 
lated to DMC in the sense that they are based on a 
drift-diffusion process in imaginary time. Unlike DMC, 
however, these methods are based on a true Metropolis 
algorithm without the need for branching. Population 
correlation does not occur and the statistical weight fluc- 
tuations are not a problem. Instead, an analysis of the 
statistical efficiency of these methods would need to take 
into account the correlation time and its scaling with sys- 
tem size. 



X. CONCLUSIONS 

To conclude, we have derived an expression for the 
scaling behavior of DMC calculations when aiming at 
a fixed statistical precision per particle. Using a linear 
scaling algorithm for an individual time step, constant 
scaling of the total computational cost for the energy per 
particle is possible in principle, except for a factor Xpop, 
which quantifies the correlation within the population of 
walkers. The exact value of Xpop was derived in the per- 
turbative limit, depending only on the correlation time 
and the raw variance of the DMC process. Based on nu- 
merical evidence, we demonstrated that an exponential 
extrapolation of the pcrturbative law gives a strict lower 
bound to the inefficiency factor Xpop- From this, it fol- 
lows that the DMC algorithm generally scales at least 
exponentially in the square root of the system size. 

The numbers for actual sample systems indicate that 
this exponential scaling should not even be observable in 
most DMC based studies done so far, leaving plenty of 
room to do interesting research with the DMC method. 

Alternative schemes for branching and population con- 
trol that have been suggeste d 23 ' 24 i 27 may certainly influ- 
ence the efficiency of the algorithm. The exponential 
lower bound of the statistical inefficiency, however, may 
at best be shifted over towards an exponentially scaling 
population control bias^ 

It must be stressed that this exponential scaling factor 
is specific to the DMC method and does not occur in 
other methods like VMC. It is not linked to the more 
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fundamental fcrmion sign problem^ and it is not limited 
to certain observables^ 

In fact, the exponential scaling may be a symptom of 
the very nature of the DMC process. In general, Markov- 
chain MC methods such as VMC exhibit excellent scal- 
ing behavior. DMC however, is not based on a Markov 
process but rather on the simulation of a time depen- 
dent stochastic diffusion process. As such, it must be 
expected to suffer from the exponential accumulation of 
errors inherent in the simulation of time evolution in non- 
integrable systems. The population control that is neces- 
sary to stabilize the process might then necessarily lead 
to exponential scaling either in the bias or the efficiency 
of the process. The only way to overcome this problem 



might then be to resort to alternative QMC methods like 
VMC, PIMC or RMC that arc based the stochastic com- 
putation of a multi-dimensional integral in the original 
spirit of Markov-chain MC methods. 
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